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Abstract. The formation of singularities in solutions to the dispersionlcss Kadomtsev-Petviashvili 
(dKP) equation is studied numerically for different classes of initial data. The asymptotic behavior 
of the Fourier coefficients is used to quantitatively identify the critical time and location and the 

ff\ type of the singularity. The approach is first tested in detail in 1 + 1 dimensions for the known case 

of the Hopf equation, where it is shown that the break-up of the solution can be identified with 

(•"■■v prescribed accuracy. For dissipative regularizations of this shock formation as the Burgers' equation 



(N 



and for dispersive regularizations as the Korteweg-de Vries equation, the Fourier coefficients indicate 
as expected global regularity of the solutions. The Kadomtsev-Petviashvili (KP) equation can be 
seen as a dispersive regularization of the dKP equation. The behavior of KP solutions for small 
dispersion parameter e <C 1 near a break-up of corresponding dKP solutions is studied. It is found 
^T that the difference between KP and dKP solutions for the same initial data at the critical point 

scales roughly as e 2 ' 7 as for the Korteweg-de Vries equation. 
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1. Introduction. Nonlinear evolution equations without dispersion and dissi- 
pation generically have solutions which show the wave breaking phenomenon, i.e., the 
formation of a shock, a gradient catastrophe in finite time. A standard example in 
i ^ i this context is the Hopf equation 



> 



u t + Quu x = 0. (1-1) 



It is well known that the solution of ( 1.1 ) for an initial value problem u(x, 0) = uq(x) 



can be obtained via the method of characteristics in the implicit form 



Jg u(M)=«o(0. a> = 6tuo(0+f- (1-2) 

If the initial data are such that t c = -. is positive, the solution reaches 

max feK [-6u[,(q] 
a point of gradient catastrophe x c at t c where the derivative of the Hopf solution blows 
up, but where the solution stays finite. 

^ Regularizations of this equation with small dissipation e, the Burgers' equation, 

• i— i 

^ u t + Quu x = eu xx , (1.3) 

H 

Cu or small dispersion e, the Korteweg-de Vries (KdV) equation, 

u t + 6uu x + e 2 u xxx — 0, (1.4) 
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will have solutions which stay regular at the shock of the Hopf solution for the same 
initial data, but show some critical behavior in the vicinity of the point (x c , t c ). At the 
critical point, the difference between Hopf solution and the solution of the regularized 
equation shows a characteristic scaling in e, for KdV e 2 ' 7 . Dubrovin [3] conjectured 
a universal behavior of solutions to Hamiltonian regularizations of the Hopf equation 
(among which KdV is the most prominent example) in the vicinity of a shock. In 
[T51 [TOl ITT1 |2"U] strong numerical evidence for this conjecture was given, which was 
proven for KdV in [5] via Riemann-Hilbert techniques. For t -C t c this difference 
scales as e 2 . An asymptotic description of dissipative regularizations was presented 
in [TO]. 

Dubrovin's conjecture is based on a double scaling limit where e — > and simul- 
taneously x — > x c , t — ¥ t c in such a way that the limits 



lim 
6u c t — > x c — 6u c t c 



x c - 6u c (t - t c ) 



and 



lim 

t ->tr 



(t - t c ) 



where u c — u(x c ,t c ), exist and are bounded. In [TH1 UHl HU I2S] it was shown that it 
is indeed possible to study these scalings in e numerically if the critical point (x c ,t c ) 
is known. The above formulae make it clear that the error in the values x c , t c must 
be smaller than the considered values of e if the correct scalings are to be identified 
numerically. This implies that for a numerical study of critical scaling phenomena, a 
control of the error bounds for the critical point is crucial. 

It is mathematically and physically interesting to study similar problems in higher 
dimensions where much less about shock formation and dissipative or dispersive reg- 
ularizations is known. A breakdown of regularity always indicates a limit of the 
applicability of the studied model. It is at such points that effects as dissipation and 
dispersion, which have been neglected in the simplified model, become important. The 
main technical problem in higher dimensions is that the dispersionless system even of 
completely integrable equations as the Kadomtsev-Petviashvili (KP) equation [21], a 
2 + 1-dimensional variant of KdV, is not integrable in the usual sense. This means 
that there are in general no standard solution generating techniques as hodograph 
methods or (linear) Riemann-Hilbert problems available in this context (but see [25] 
for a nonlinear Riemann-Hilbert approach). The KP equations read 



d x (d t u + 6ud x u + e 2 d xxx u) + \d yy u = 0, A = ±1 



(1.5) 



where (x,y,t) £ Rj x 1 9 x I ( and where e <C 1 is a small scaling parameter. The 
case A = — 1 corresponds to the KP I model with a focusing effect, and the case A = 1 
corresponds to the KP II model with a defocusing effect. 

The focus in this paper will be on the dispersionless variant of the KP equation, 
the dKP equation, also known as Khokhlov-Zabolotskaya equation [46 , which follows 



from (1.5| for e = 0, 



d x (d t u + 6ud x u) + XdyyU = 0, A = ±1. 



(1.6) 



It appears in many applications as a model for dissipation- and dispersionless, essen- 
tially one-dimensional waves with weak transverse effects, for instance in nonlinear 
acoustic and in gaz dynamics, see [35J[5T]. It also plays a role in the context of general 
relativity and differential geometry [TJ [T2]. Alinhac and coworkers, see for instance [5] 
and references therein, showed that it appears as a universal model in the geometric 
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study of singularity formation in nonlinear wave equations. Existence of solutions up 
to break-up was proven, and the singularity was identified for a generic setting as a 
one-dimensional cusp. 

The dKP equation is not completely intcgrablc in the sense that an infinite number 
of conserved quantities exists (in fact it only has three) . But it is integrable according 
to the definition of [13] that an infinite number of hydrodynamical reductions exists, 
see also [26, 271133 • As was shown in [T^], solutions can be constructed with Twistor 
methods in terms of Einstein- Weyl metrics. In [30] EH] solutions to the dKP equation 
were obtained via the solution of a nonlinear Riemann-Hilbert problem. This allowed 
the study of the long-time asymptotics and finite time break-up in J3TJ [25J . Another 
way to construct dKP solutions arises from the theory of Frobenius manifolds. In [57] 
the infinite dimensional Frobenius manifold corresponding to dKP was constructed. 
A common feature of the above solution generating techniques is that even the most 
explicit forms of the solutions require the solution of singular nonlinear integral equa- 
tions and the inversion of implicit functions which makes it difficult to find explicit 
examples for the wave breaking. Since none of these approaches has been studied 
numerically so far, we solve here dKP directly up to the formation of a singularity. 

As mentioned above, a precondition for the numerical study of scaling laws as 
for KdV is that the critical point and the critical solution can be determined with 
sufficient accuracy. It is the goal of this paper to provide the necessary tools for 
this in a rather general setting, and to apply them to the dKP equation for several 
classes of initial data. The applicability of this approach will be shown first for the 
1 + 1-dimensional case where exact solutions provide tests. Note that a dissipative 
regularization does not give precise information on the critical point as can be inferred 
from the example of the Burgers' equation with small dissipation we will also consider. 
The same will be shown to be true for a dispersive regularization as provided by KdV. 
Thus it can be concluded that a direct study of the dispersion- and dissipationless 
equation with the techniques explored in this paper is more promising in the context 
of identifying singularity formation. 

The basic idea of the used numerical approach is to approximate the spatial de- 
pendence of the solutions by a discrete Fourier series. It is known that the asymptotic 
behavior (for large wave numbers k) of the Fourier coefficients is exponential for solu- 
tions analytic in the complex plane in the vicinity of the real axis. This dependence 
becomes algebraic at the break-up of the solution. Thus the vanishing of the ex- 
ponential decrease in the Fourier coefficients would indicate the appearance of the 
singularity. In practice there are several problems in this context: firstly it is difficult 
numerically to integrate the equation up to the critical point; moreover even if this can 
be done with sufficient accuracy, the Fourier coefficients might be polluted especially 
at the high wave numbers where the asymptotic behavior is to be read off; and last 
but not least, the asymptotic relations for the Fourier coefficients are established for a 
continuous Fourier transform, whereas numerically only a discrete Fourier transform 
is considered as an approximation. 

Asymptotic Fourier analysis was first applied numerically by Sulem, Sulem and 
Frisch [31] , at the time of course with much lower resolution than is accessible today. 
Therefore the study in [41] was mainly qualitative, but gave a proof of concept. In 
[42] they studied singular solutions to the two-dimensional cubic NLS equation. An 
application of the method to the 2-d Euler equations can be found in [T7J [35] . It has 
also been applied to the study of complex singularities of the 3-d Euler equations in 
[3], in thin jets with surface tension [3j5], the complex Burgers' equation [32] and the 
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Camassa-Holm equation [38j to name just a few examples. The main goal in these 
works is to decide whether the studied solutions develop a singularity. However, the 
purpose of the present paper is to identify time and location of a singularity, which 
is known to appear in the solution, with prescribed accuracy with these methods. 

The tracking of singularities can be of course also achieved with other techniques, 
for instance Pade approximants in [45]. The problem with the latter is, however, 
that this adds an additional dimension to the problem which is especially expensive 
if this technique is to be applied in higher dimensions. Full singularity tracking for a 
function in K™ would imply the study of the function in C™, thus doubling the spatial 
dimensions. Therefore a quantitative numerical approach based on the asymptotic 
behavior of discrete Fourier series would be an economic way to study singularity 
formation in higher dimensions for functions analytic before a critical time t c . It is 
one goal of this paper to present a careful investigation of the possibilities and the 
limitations of this approach, and to show that the method can provide information 
on critical points with the needed accuracy to study scaling laws. The techniques are 
developped for the example of the Hopf equation and various regularizations thereof, 
and will then be applied to the dKP equation for various initial data. The thus 
identified dKP solution at the critical point is used to study the scaling in the small 
dispersion parameter e of corresponding KP solutions. It is found that the scaling of 
the difference between KP and dKP solutions is compatible with e 2 / 7 . This suggests 
that the behaviour of the KP solutions in the singular direction is similar to KdV 
solutions near the break-up of Hopf solutions. 

To establish asymptotic Fourier analysis as a quantitative tool to identify shock 
formation, we test the various numerical problems in the approach |41) separately and 
show they can be controlled. Then we will establish that this method in fact can be 
used to determine the critical point and the critical solution to the Hopf equation. 
The paper is organized as follows: In Sec. 2 we collect some known facts about the 
asymptotic behavior of Fourier coefficients of a real function analytic in the vicinity of 
the real axis. We numerically integrate the Hopf equation up to the critical point for 
given initial data and compare it with the exact solution. For the latter we show how 
well the Fourier coefficients can be fitted to the expected asymptotic decrease. We 
perform this fitting during the numerical solution of the Hopf equation to determine 
the critical time. The same approach is applied to solutions to the Burgers's and 
KdV equation for the same initial data. In Sec. 3 we apply these techniques to the 
break-up in dKP solutions for two classes of initial data, localized in two dimensions 
or as the line solitons of KP, localized in one spatial direction and infinitely extended 
in the other. In both cases, the scaling of the difference between KP solutions in the 
small dispersion limit and dKP solutions at the critical point is studied. We add some 
concluding remarks in Sec. 4. 

2. 1+1-dimensional case. In this section we will address numerically 1 + 1- 
dimensional examples from the family of Hopf equations. We first review some well 
known facts from asymptotic Fourier analysis. Then we will numerically integrate 
the Hopf equation for a concrete example up to shock formation. We discuss how the 
asymptotic behaviour of the Fourier coefficients can be used to identify the critical 
time and solution of the Hopf equation. A similar analysis is then performed for the 
Burgers' and the Korteweg-de Vries equation which can be seen as dissipative and 
respectively dispersive regularizations of the Hopf equation. 

2.1. Asymptotic Fourier analysis. In this subsection we review some well 
known facts on the asymptotic behavior of the Fourier transform of a real function 
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which is analytic in a strip around the real axis. The resulting formulae will be then 
used to numerically track singularities of the solutions in the complex plane via their 
Fourier coefficients without having to deal with two real dimensions as in |45) . This 
allows to determine when a singularity hits the real axis, i.e., when the real solution 
becomes singular. 

Spectral methods are a powerful tool in the numerical solution of differential 
equations because of their excellent approximation properties for smooth functions, 
see for instance [5]. It is well known that the error in the approximation of such 
a function with a spectral approximant of order N, e.g. a polynomial of degree iV, 
decreases faster than any power of 1/N. This provides a huge advantage over methods 
with just an algebraic decrease of the error and is the basis of the efficiency of spectral 
methods, see also [IS]- The most used such methods are the expansion of the studied 
function in terms of discrete Fourier series. For given analyticity properties of the 
studied function, precise mathematical statements on the asymptotic behavior of the 
Fourier coefficients exist. As we will detail below, these imply that the location 
of singularities in the complex plane can be obtained from a given Fourier series 
computed on the real axis. It should be possible in this way to numerically determine 
the type of the singularity, to trace it in the complex plane, and to establish when 
the singularity hits the real axis as was done for the first time in [4"T] . 

The Fourier transform u of u(x) <G L 2 (R) is defined as 

u{k) = f u{x)e~ lkx dx. (2.1) 

A singularity in the complex plane of the form u ~ (z—Zj)^ , \ij <£ Z, with Zj = aj—i5j 
in the lower half plane (8j > 0) implies with a steepest descent argument for k — > oo 
the following asymptotic behavior of the Fourier coefficients (for a detailed derivation 
see e.g. [5]), 






2n^~ 2 e-^ y ' g-tucti-M^ (2.2) 



Consequently for a single such singularity with positive dj , the modulus of the Fourier 
coefficients decreases exponentially for large k. For Sj — 0, i.e., a singularity on the 
real axis, the modulus of the Fourier coefficients has an algebraic dependence on k. If 
there are several singularities of this form at Zj, j = 1, . . . , J, there will be oscillations 
in the modulus of the Fourier coefficients for moderately large k. 

To numerically compute a Fourier transform, it has to be approximated by a 
discrete Fourier series which can be done efficiently via a fast Fourier transform (FFT), 
see e.g. |43j . The discrete Fourier transform of the vector u with components Uj — 
u(xj), where Xj = 2irLj /N , j = 1, . . . , N (i.e., the Fourier transform on the interval 
[0, 2irL] where L is a positive real number) will be always denoted by v in the following. 



There is no obvious analogue of relation (2.2) for a discrete Fourier series, but it can 



be seen as an approximation of the latter, which is also the basis of the numerical 
approach in the solution of the PDE. It is possible to establish bounds for the discrete 
series, see for instance [3J. We will in this paper always fit the discrete Fourier series 



for large \k\ to the asymptotic formula (2.2) 



Remark 2.1. Note that the minimal distance in Fourier space is m :— 2ttL/N 
with N being the number of Fourier modes and 2-kL the length of the computational 
domain in physical space. This defines the smallest distance which can be resolved in 
Fourier space. All values of 5 below this threshold cannot be distinguished numerically 
from 0. 
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2.2. Numerical integration of the Hopf equation. In this subsection we 
will study the integration of the Hopf equation for a concrete example up to the 
appearance of a gradient catastrophe. Since the method of characteristics gives an 
exact solution in this case in implicit form which can be obtained numerically in 
principle with machine precision, we can use this as a test case for the numerical 
integration of a nonlinear dispersionless equation up to the first critical point in the 
solution. We discuss how to obtain asymptotically reliable Fourier coefficients, and 
with which precision this can be achieved. 

In the following we will consider the example uq(x) = sech x for which the gra- 
dient catastrophe appears at the critical time t c = \/3/8 ~ 0.2165 and at the critical 
point x c = v3/2 — ln((v3 — l)/v2) ~ 1.5245. These initial data are motivated 
by the well known KdV soliton. They have the advantage that they belong to the 
class of rapidly decreasing functions which, for sufficiently large L, can be analytically 
continued as periodic functions within numerical precision (the quantity L and thus 
the length of the interval [— ttL,ttL] is chosen large enough that uq(±itL) is smaller 
than machine precision (10 -16 in our case). This ensures that there is no Gibbs phe- 
nomenon and that the analytic behavior of the Fourier coefficients is not affected by 
effects at the boundary of the computational domain. 



Numerically we solve equation (1.2 1 iteratively 



e 



n+l 



6tUo(£ n 



with the initial guess £o = x. If some relaxation is used, the iteration converges and 
is stopped once the residual of 6tu (£) + £ — x is smaller than 10~ 13 in the L^ norm. 
The iteration is done in a vectorized and thus efficient way, i.e., for all points at the 
same time. For the studied example the solution to the Hopf equation can be seen in 
Fig. 



2.1 at the initial and at the critical time. 




Fig. 2.1. Solution to the Hopf equation (l.ly propagating to the right for uo(x) 
t = in blue and at the critical time t c = v3/8 ~ 0.2165 in green. 



sech x at 



To directly integrate the Hopf equation numerically, we use a standard Fourier 
spectral method in x: for x G [— 7r,7r]L, we take L = 5, and the modulus of the 
Fourier coefficients computed via a discrete Fourier transformation decreases to ma- 
chine precision for more than N = 2 8 Fourier modes. Treating the spatial dependence 



in (1.1) with a discrete Fourier transform, we obtain for the Hopf equation a system 
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of ordinary differential equations (ODEs) in t for the Fourier coefficients. Since the 
latter system is not stiff in contrast to the KdV equation treated in a later subsection, 
we can use standard time integrators here. We choose the well-known explicit fourth 
order Runge-Kutta scheme for convenience. 

The accuracy of the numerical solution is tested via comparison with the exact 
solution. In addition we trace a conserved quantity of the Hopf equation, the energy, 



E[u] = / u 3 dx, 



which, when computed for the numerical solution, will depend on time due to unavoid- 
able numerical errors. As shown for instance in [231 IH], the conservation of energy 
can be used as an indicator of numerical accuracy. We test this for the Hopf equation 
below by tracing the quantity AE = |1 — g ["J,Q, |, where E[u](t) is the numerically 
computed energy. 

Since the solution of the Hopf equation is very different for times t <C t c and for 
t ~ t c , we will first numerically solve the Hopf equation up to t = 0.18 <C t c . For 
N = 2 14 Fourier modes, the solution at time t = 0.18 is fully resolved as can be seen 



from Fig. |2.2| where the modulus of the Fourier coefficients is shown. The modulus 
goes down to machine precision. For a time step A t = 7.2 * 10~ 5 the L^ norm of the 
difference between numerical solution and exact solution is 6.4 * 10 , the quantity 
AE ~ 10~ 14 . For At = 3.6*10 -5 , the respective numbers are ||u„ um — u exac t\\oo = 3.7* 
10~ 12 and AE ~ 10~ 15 , for A t = 1.8 *10" 5 one finds ||u„„ m -u ex0 ctlU = 8.4*10~ 13 
and AE ~ 10~ 14 . Thus AE overestimates the numerical precision by roughly two to 
three orders of magnitude (once it is of the order of machine precision, it is obviously 
dependent on rounding errors and can even increase with smaller resolution as above). 
It can be in fact used as an indicator of numerical accuracy in cases for which no exact 
solution is known, provided that there is sufficient spatial resolution. To sum up, the 
solution can be obtained up to machine precision for times t <^ t c . 




-2000 -1500 -1000 -500 



Fig. 2.2. Modulus of the Fourier coefficients for the solution to the Hopf equation \1.1\ fo 
uq(x) = sech?x at t = 0.18. 



For larger t, time steps of different size might be convenient. Thus we run the 
code with a certain time step up to time t = 0.18 and use the result as initial data 
for a second time evolution up to the critical time t c (since we have an exact solution 
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for the studied example, it could be also used as exact initial data at that time; but 
since the goal is to develop and test methods for a general case for which no exact 
solution is known, this is not done here). The important point is that the solution 
can be obtained with machine precision at any time tct c numerically. 

For the studied example, one gets for N — 2 14 with A t = 7.3 * 10 -5 an L^ norm 
of the difference between numerical and exact solution of the order of 0.032. This 
error (as well as the L2 norm of the difference between numerical and exact solution) 
does not decrease for smaller time steps. This shows that the error, which is always 
biggest at the critical point as can be seen in Fig. |2.3| is due to a lack of resolution in 
x and not in time. This lack of resolution is clear from the Fourier coefficients which 



can be also seen in Fig. 2.3 Visibly the Fourier coefficients for high wave numbers do 
not show the expected decrease with k due to numerical errors at the critical time. 
This will be important for the asymptotic analysis in the following section. 
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Fig. 2.3. Difference of the numerical and exact solution at time t c for N = 2 14 Fourier modes 
on the left, and the modulus of the Fourier coefficients for the numerical solution on the right. 



Note that the quantity AE for the computed energy decreases to machine preci- 
sion if the time step decreases to A t — 9.13 * 10 -6 . This shows that this quantity is 
only a useful indicator of the numerical accuracy for sufficient spatial solution, i.e., if 
the Fourier coefficients decrease to the wanted accuracy. Lower values of AE than the 
modulus of the Fourier coefficient of the highest wave numvers are meaningless. To 
obtain higher accuracies, it appears necessary to allow for higher spatial resolution. 
For N = 2 15 Fourier modes and A t = 3.7 * 10~ 5 , we get ||u„ um — u exact || 00 = 0.025, 
for N = 2 16 and A t = 1.8 * 10~ 5 ||u rmm — u e:m ct||oo = we find 0.016. The time steps 
have been chosen in a way that the error does not decrease further if A t is decreased. 

2.3. Asymptotic fitting of the Fourier coefficients. In this subsection we 
will study the asymptotic behavior of the Fourier coefficients for the solution to the 
Hopf equation discussed in the previous subsection. We show how to apply the asymp- 
totic formula (2.2 1 to the exact solution before and at the critical time. This will be 



done here first for the exact solution to avoid the problems with errors in the numer- 
ically determined Fourier coefficients at the critical time for high wave numbers. 



As in the previous subsection, we discuss the Hopf equation (1.1) for the initial 



data Uq(x) — sech x. For a time t -C t c , the solution and the Fourier coefficients 
computed via FFT with N = 2 14 Fourier modes on the interval x £ [—ir,ir]L with 
L = 5 can be seen in Fig. |2.4| The latter show the expected exponential decrease up 
to a level of roughly 10~ 13 where the error saturates because of the finite numerical 



precision. 
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Fig. 2.4. Exact solution of the Hopf equation for the initial data uq(x) 
t c on the left, and the modulus of the Fourier coefficients on the right. 



sech 2 x fort = 0.18 < 



In this case, only one singularity is expected in the lower half plane. To test the 



asymptotic formula (2.2) in practice, we first do a least square fitting for ln|v|, which 
should be according to (2.2) of the form 



]n\v\ ~ A- Bhik- k5. 



(2.3) 



The fitting is done for a given range of wave numbers k (we only consider positive 
k). For obvious reasons, k has to be limited to values for which \v\ is larger tha n the 
rounding errors. We choose this threshold to be |«| > 10 -12 . In principle formula (2.2 ) 
holds only for k S> 1. However if we do the fitting for all positive k, we get A — 6.6874, 
B = 1.4309 and 5 — 0.0349 as well as A = 2.05, where we define the quantity A as 
the Lao norm of the difference between the solution and the the asymptotic formula 



(2.3) 



hi\v\ - (A -B hi k-k8)\\ 



(2.4) 



As expected the difference is maximal near k = 0. For k m i n = 10 we get A — 6.68693, 
B = 1.4731 and 8 = 0.0348 and 0.0293 for the quantity A. The difference between 
In \v\ and the fitted formula is shown in Fig. 2.5 The main difference appears for small 
k. From the figure it can be concluded that a minimal value of k m i n = 100 might be 
optimal. On the other hand it can be also seen that rounding errors lead to some fuzzy 
behavior of the Fourier coefficients at high wave numbers. This suggests to consider 
only values \v\ > 10~ 10 . With these choices, we obtain A = 6.9311, B — 1.4863 
and S = 0.0347. The difference between In \v\ and the fitted formula can be seen in 
Fig. 2.5 it is of the order 10 -4 . Note that the value of S is not much affected by the 
different choices of these thresholds. This just reflects the fact that the fitting to an 
exponential is much more stable than the fitting to an algebraic function. 

The Loo norm A (2.4 ) of the difference between In \v\ and the fitted formula (2.2 ) 



can be used to indicate the quality of the fitting. An additional test is provided 
from the relation between A and B. Since we compute the Fourier coefficients via 



a discrete Fourier transform whereas formula (|2.2|) holds for a standard continuous 
Foi 

(B 



Fourier transform, an additional term In ^^ appears in the formula, A = \ ln(27r) 



ln(B - 1) - (B - 1) + In 



2-kL 

N ■ 



We get for the right hand side of this relation 
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Fig. 2.5. Difference between t he m odulus of the Fourier coefficients for the situation of Fig. \2Ji\ 
and the fitted asymptotic formula \2.efy for k > 10 and \v\ > 10~ 12 on the left, and for k > 100 and 
\v\ > 10~ 10 on the right. 



6.5352, a value close to the computed A — 6.9311. However it is clear that the 
quantity A is the more convenient indicator of the quality of the fitting. Thus it will 
be used to this end in the following. 

The above fitting has only been carried out for the absolute value of the Fourier 
coefficients. To determine a, the real part of the location of the singularity, we fit the 
imaginary part of the logarithm of v, 



:= 3dni> ~ C — ak. 



(2.5) 



Since the logarithm is branched in Matlab at the negative real axis with jumps of 27r, 
the computed <f> will in general have many jumps. Thus one has first to construct 
a continuous function from the computed <p. The analytic continuation will be done 
in the following way: starting from the first value (largest fc), we check for all other 



values oi(j)(kj) whether \4>(kj +1 ) — 4>{kj)\ > \cj){k 



j + i) — (f)(kj)±ir\. If this is the case, we 
put cj)(kj + i) — > </>(fcj+i)±7r. The result of this procedure will be a continuous function 
which will be fitted with a least square approach to a linear function. For the above 
example, we get a = 1.3754 and C = 136.5425. The L^ norm of the difference of the 
function <f> and the fitted straight line, 

A 2 = \\%\nv-C + ak\\ 00 , 

can be seen in Fig. |2.6| It is of the same order of magnitude as the quantity A which 
shows the self consistency of the approach. 

The above situation is typical for a single singularity in the complex plane away 
from the real axis. This picture changes if S = 0, i.e., when the singularity hits the 
real axis which happens for the considered example at the critical time t c . In this case 
the Fourier coefficients as expected no longer show an exponential decrease as can be 
seen in Fig. |2.7| Note the difference to the numerically computed Fourier coefficients 



at high wave numbers in Fig. 2.3 



The asymptotic fitting of \n\v\ as in Fig. [2_5j for k > 100 yields A = 7.0283, 
B = 1.4249 and 5 = —0.0001. The difference between In \v\ and the fitting curve can 
be seen in Fig. |2.8| Since there is no exponential decay in this case, the fitting is 
much more sensitive to the considered cutoffs, for low k due to the condition k ^$> 1 
for the asymptotic formula to be valid, for large k due to rounding errors. This is 
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Fig. 2.6. Differenc e Ag between the imaginary part of the loga rithm of the Fourier coefficients 
for the situation of Fig. \2.J\ and the fitted asymptotic formula \2.2\ for k > 100 and \v\ > 10~ 10 . 
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Fig. 2.7. Modulus of the Fourier coefficients of the exact solution of the Hopf equation for the 
initial data uq(x) = sech x at the critical time t c . 



also obvious from Fig. |2.8| where the errors are maximal at both ends of the shown 
spectrum. For 200 < k < 900 we get A = 6.4573, B = 1.3084 and S = 0.0001 with a 
fitting error as shown in Fig. 2.8 In both cases 5 is of the order of 10~ 4 . The quantity 
B should be equal to 4/3 and is as before much more sensitive to a choice of the limits 
kmin and kmax for the fitting. 

By fitting the quantity <f> for 200 < k < 900, we get a = 1.5251 and C = 306.3523 
with a fitting error of the same order. The former is very close to the exact value 
x c = ^3/2 - ln((v / 3 - l)/-\/2) ~ 1.5245. If we repeat the above fitting for N = 2 15 
Fourier modes for 200 < k < N/3/L (L = 5 as before), we obtain A = 7.4012, B = 
1.3597 and 5 = -3.9054 * 10~ 5 with a A = 0.012 and a = 1.5248 and C = 306.3584. 
Thus a higher number of Fourier modes makes the fitting more stable and reliable. 

To sum up, we have shown in this subsection that the quality of the fitting depends 
on the choice of the minimal and maximal values for the wave numbers. Generally 
a choice of k m i„ = 100 or higher appears to be recommended. However this has the 
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Fig. 2.8. Difference between t he m odulus of the Fourier coefficients for the situation of Fig.iK 
and the fitted asymptotic formula tfTw for k > 100 on the left, and for 200 < k < 900 on the right 



disadvantage that many Fourier modes are excluded, especially at times i<t c where 
only few modes contribute. Rounding errors enforce a choice of k max of the order of 
N/3/L. By studying the difference between the considered quantity and the fitting 
curve, these values can be optimized as outlined in the next subsection. A proper 
choice of these thresholds is more important at the critical time when the singularity 
hits the real axis since there is no exponential decay in this case. 

2.4. Asymptotic behaviour of the Fourier coefficients for a numerical 
solution to the Hopf equation. In the previous subsection we have studied the 
asymptotic fitting of the Fourier coefficients on a finite interval of the wave numbers 
and have fixed the limits of this interval essentially by hand. Here we will do this 
in a more systematic way by choosing a suitable lower threshold (which depends on 
the studied problem) , and by varying the upper limit to reach a prescribed difference 
between the Fourier coefficients and the fitted curve with a maximal number of co- 
efficients. Since this difference would obviously be minimal if upper and lower limit 
coincide, we impose the additional condition that at least half of the Fourier coeffi- 
cients with \v\ > 10 -10 will be included in the fitting to include most of the available 
information on the solution. With this approach and the preliminary studies of the 
previous sections, we are then able to perform the asymptotic fitting during the ac- 
tual numerical computation of the Hopf solution. We show that the vanishing of the 
computed 5 can be in fact used to determine the critical time t c and the type of the 
singularity via the exponent fj,. The goal is to obtain an error in the fitting of the 
same order as the numerical error. 

The fitting method used in this paper is obviously not the only possible one. 
For instance in [4] and related publications the sliding fit approach is applied where 
only few values of k are taken into account in the procedure (typically between 4 
and 100 for a resolution of 1024 modes, see for example [35] )■ The quality of the 
fitting in this approach is roughly determined by the independence of the computed 
fitting parameters on the minimal value k m i n of the k used in the fitting, and on the 
discretization size 2ttL/N. Note that the sliding procedure confirms the results found 
here, for k m i n < 100 and up to 600 points considered. Nevertheless, this approach is 
computationally expensive which implies that the time of the singularity formation 
has to be previously roughly determined in another way to be able to perform the 
study close to it. Since it is the aim of this paper to identify the break-up numerically 
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via the asymptotic behavior of the Fourier coefficients, even in cases where the critical 
time is not known, we do not use the sliding fit here. 

As in the previous section, the asymptotic fitting of hi(|v|) is done on the interval 
k m in < k < k max , where we put 



10 and k.„ 



min^!,^), with 



k\ = max (k), ki — max(fc) 

|u|>le— 10 A<p 



where p denotes a prescribed value for A in (2.4) on the fit interval. The choice 
kmin — 10 allows the fitting even for short times where there are only few Fourier 
coefficients with \v\ > 10~ 10 . It satisfies the condition k m i n 3> 1 whilst maximizing 
the number of the numerically computed Fourier coefficients included in the fitting 
procedure. Indeed, for given p, larger values of k m i n would lead to a too small number 
of coefficients used for the fitting. 

For N = 2 15 , we observe the influence of p on the time of vanishing of 5 and the 
corresponding value of B in Fig. |2.9| and in Table [2~T| Due to the large computational 
domain considered (x G [— w,iv]L, L = 5), a sufficient resolution has to be used to get 
a sufficiently small value of m (m := 2irL/N, m ~ 10~ 3 for N = 2 15 and L = 5). 




0.21 0.211 0.212 0.213 0.214 0.215 0.216 



Fig. 2.9. Time evolution of the fitting parameters S and B for the numerical solution to the 
Hopf equation with u(x, 0) = sech 2 (x), with N = 2 15 and different values of p between t = 0.21 and 
t c = V3/8. The dashed line in the left figure corresponds to the minimal resolution in Fourier space 
discussed in remark] 2. 1\ 



But for smaller 



If p decreases to 10 3 , one finds that (t c ,B(t c )) — >• ( g , 

values of p, for example for p — 5 * 10~ 4 , the values for both t c and a(t c ) become 
worse approximations of the exact ones. In fact, to get such small fitting errors, we 
have to consider too few Fourier coefficients which explains the lack of precision for 
these. Thus we add the requirement that at least half of the total number of the 
Fourier coefficients available have to be used for the fitting. For N = 2 15 , this implies 
that the 'minimal' achievable fitting error whilst using enough Fourier coefficients is 
p = 0.01, which corresponds roughly to the precision with which we compute the 
numerical solution, see Sec. |2.2| 

Thus for different values of N, and the related precision with which we compute 
the numerical solution, we observe in Fig. 2.10 the time evolution of 8 and B, and 

As N increases, 
(tc,B(t c )) 



give the values of 5, B and a at two different times in Table 2.2 

( ^, | ) and a(t c ) — > 1.5245 without further restrictions on k. The time 
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P=1Q- / 


t 


5 


B 


a 


£<5<0 


0.2136 
0.2164 


9*10~ 4 
-5*10~ 6 


1.363 
1.339 


1.5128 
1.5242 




p = 5 * 10" 3 


t 


S 


B 


a 


ts<o 


0.2138 
0.2165 


9*10~ 4 
-7*10~ 6 


1.356 
1.337 


1.5134 
1.5244 




p=10- 6 


t 


8 


B 


a 


ts<o 


0.2143 
0.2165 


9*10~ 4 
-7*10- 6 


1.345 
1.334 


1.5149 
1.5245 




p = 5*10~ 4 


t 


S 


B 


a 


ts<o 


0.2145 
0.2166 


9*10~ 4 
-1*10" 5 


1.342 
1.333 


1.5153 
1.5243 



Table 2.1 
Values of the fitting parameters for the numerical solution to the Hopf equation with u(x, 0) = 
sech(x) 2 , with N = 2 15 and different values of p = 0.01,0.005,0.001,0.0005. The exact values are 
t c ~ 0.2165, B = 4/3 and x c ~ 1.5245. 




Fig. 2.10. Time evolution of the fitting parameters 8 and B for the numerical solution to the 
Hopf equation with u(x,0) = sech(x) 2 , for different values of N between t = 0.21 and t c = \/3/8. 



dependence of the fitting parameters is very similar for different N, which indicates 
that the fitting is reliable. 

We conclude from this example that the method can be used in practice to determine 
the critical time t c and the type of the singularity via the value of B(t c ) if sufficient 
spatial resolution is used. The choice of the minimal wave number k m i n for the fitting, 
and the achievable minimal fitting error depend on the studied problem. 



2.5. Burgers' equation. The Burgers' equation (1.3 1 with small dissipation 



< t « 1, can be viewed as a purely dissipative regularization of the Hopf equation. 
An asymptotic description at the critical point was given in |10j . Since the system of 
ODEs resulting from the approximation of the solution via a discrete Fourier series in 
the spatial coordinate is now stiff, we use a stiff integrator, a fourth order exponential 
time differencing scheme by Cox and Matthews [7] as implemented in 39, 23 . This 
method is very efficient for both Burgers' and KdV equation, see e.g. 23 . In both 



cases we consider as before the initial data Uq{x) 
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sech x for x e [— 5n, 57rl 



N=2 U 


t 


S 


B 


a 


£<5<0 


0.2116 
0.2164 


0.0019 
-6*10" 6 


1.37 

1.345 


1.5048 
1.5238 



AT = 2 15 


t 


6 


B 


a 


£<5<0 


0.2136 
0.2164 


9 * 10- 4 
-5*10" 6 


1.363 
1.339 


1.5128 
1.5242 



iV = 2 16 


t 


S 


B 


a 


£5<m 

ts<o 


0.2148 
0.2165 


4.7*10-4 
-5*10" 6 


1.355 
1.338 


1.5176 
1.5244 



iV = 2 17 


t 


6 


B 


a 


ts<o 


0.2155 
0.2165 


2.3*10- 4 

-7*10- 7 


1.349 
1.337 


1.5204 
1.5244 



Table 2.2 
Values of the fitting parameters for the numerical solution to the Hopf equation with u(x, 0) 
sech(x) 2 , for different values of N = 2 14 , 2 15 , 2 16 and 2 17 . 



For the Burgers' equation the shock of the Hopf solution for t > t c will be regu- 
larized as can be seen for e = 0.01 in Fig. |2.11| at time t — 0.23 > t c . There is a zone 
with a strong, but finite gradient of order 1/e. 
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Fig. 2.11. Solution to the Burgers' equation for the initial data Uo(x) = sech x fort = 0.23 > t c 
and e = 0.01 on the left, and the modulus of the corresponding Fourier coefficients on the right. 



The Fourier coefficients for this solution are also shown in Fig. |2.11| It can be 
seen that they decrease (essentially exponentially) to machine precision. Doing a least 
square fitting for the modulus of the Fourier coefficients with \v\ > lO -10 for k > 100 
(this is the result of the procedure outlined in the previous section with a fitting error 
of the order 10- 2 ), we find A = 2.4858, B = 0.0273 and 5 = 0.0183. This shows that 
the parameter 5 indicating the distance between the nearest singularity to the real 
axis stays finite even for t > t c as expected. The difference between the modulus of the 
Fourier coefficients and the fitted asymptotic formula is of the order of 10- 3 . Fitting 
the quantity <f> in (2.51, we get a = 1.5778 and C — 155.5075. This implies that the 



real part of the location of the singularity is shifted towards larger values than for the 
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critical value x c of the Hopf solution which is not surprising since the Hopf shock also 
propagates with t. The difference between <fi and the fitted asymptotic formula is of 
the order of 1CP 3 , too. 

For times t 3> t c the solution in the vicinity of the Hopf shock becomes steeper 
as shown in Fig. |2.12| which means that more Fourier modes have to be used in 
order that the modulus of the coefficients decreases to machine precision (N = 2 15 
instead of N — 2 14 before). A least square fitting as before for the modulus of the 




-4000 



Fig. 2.12. Solution to the Burgers' equation for the initial data uq(x) = sech x for t = 0.4 3> t c 
and e = 0.01 on the left, and the modulus of the corresponding Fourier coefficients on the right. 



Fourier coefficients for k > 100 gives A = 3.2397, B = 0.0265 and 5 = 0.0117. The 
difference between the modulus of the Fourier coefficients a nd t he fitted asymptotic 
formula is of the order of 10 -2 . Fitting the quantity <f> in (2.5 1, we get a — 2.1657 
and C = 218.3404 with a similar quality of the fitting. 

If we trace the quantities a and 8 during the computation as a function of time 
for k > 10 (for k > 100 there are not enough Fourier coefficients above the level 
of rounding errors for small times) , we get the expected behavior as can be seen in 
Fig. |2.13| The imaginary part 5 of the singularity decreases strongly with t before 
the critical time of the Hopf solution and comes close to the axis there. After t c the 
decrease becomes very slow, and the singularity stays at a finite distance from the 
real axis. This shows that the time t c cannot be really inferred from the Burgers' 
equation with small dissipation in a reliable way. The transition in all characteristic 
quantities of the solution is not distinguished at this time, and even for t > t c the 
quantity d continues to decrease. Therefore to identify t c and x c , it appears better to 
study the Fourier coefficients of the Hopf solution without dissipation to determine 
the critical point. The real part of the singularity a can also be seen in Fig. |2.13| It 
describes in some sense the motion of the dissipative shock. 



2.6. Korteweg-de Vries equation. The KdV equation (1.4 1 with small dis 



persion e«1 can be seen as a purely dispersive regularization of the Hopf equation. 
For the initial data studied there, there will be a zone of rapid modulated oscilla- 
tions forming near the shock of the corresponding Hopf solution for t greater than the 
critical time t c as can be seen for e = 0.01 in Fig. |2.14| at time t = 0.23 > t c . 

The Fourier coefficients for this solution are also shown in Fig. |2.14| It can be 
seen that they decrease (essentially exponentially) to machine precision. Doing a least 
square fitting as before for the modulus of the Fourier coefficients with \v\ > 10~ 10 
for k > 100 (no further restrictions), we find A = 0.7638, B = 0.8165 and 5 = 
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Fig. 2.13. Fitting parameters for the solution to the Burgers' eguation for the initial data 
uq(x) = sech x for e = 0.01 in dependence on the time for 5 on the left, and a on the right. 
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Fig. 2.14. Solution to the KdV equation for the initial data uq(x) = sech 2 x for t = 0.23 > t c 
and e = 0.01 on the left, and the modulus of the corresponding Fourier coefficients on the right. 



0.0274. This shows that the parameter 5 indicating the distance between the nearest 
singularity to the real axis stays finite even for t > t c as expected. The difference 
between the modulus of the Fourier coefficients and the fitted asymptotic formula 
is shown in Fig. |2.15| It can be seen that there are damped oscillations for small 
k which die out for large \k\ and which indicate the presence of several singularities 
in the complex plane in the vicinity of the real axis. A fitting to more than one 
singularity is problematic given the limited number of Fourier coefficients. Our study 
is here of a more qualitative nature. 

Fitting the quantity <f> in ( 2.5 ), we get a = 1.5517 and C — 153.9262. This implies 



that the real part of the location of the singularity is shifted towards larger values than 
for the critical value x c of the Hopf solution which is again not surprising since the 
Hopf shock propagates. The difference between <fi and the fitted asymptotic formula 
is also shown in Fig. 2.15 There are damped oscillations in this difference for small 
k, too. 

For times ( > t c there are many oscillations in the oscillatory zone as shown 
in Fig. 



2.16 The higher number of oscillations also implies more oscillations in the 



Fourier coefficients as can be seen in Fig. 2.16 



Consequently a least square fitting as before for the modulus of the Fourier coeffi- 
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Fig. 2.15. Difference between the mod ulus of the Fourier coefficients for the situation of 
Fig. 2.1J\ and the fitted asymptotic formula \2.2\ for k > 10 on the left, and for the quantity <f> 
12. bT) on the right. 
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Fig. 2.16. Solution to the KdV equation for the initial data uq(x) = sech 2 x for t = 0.4 2> t c 
and e = 0.01 on the left, and the modulus of the corresponding Fourier coefficients on the right. 



cients for 100 < k < 1930 (to exclude the last big oscillation in the Fourier coefficients 
which is presumably a numerical artifact), we find A = 0.7981, B = 1.5684 and 
5 = 0.0164. The difference between the modulus of the Fourier coefficients and the 
fitted asymptotic formula is shown in Fig. |2.17| In contrast to the case t ~ t c shown 
in Fig. |2. 15] where there are essentially just two big oscillations, there are many more 
here. For the modulus of the Fourier coefficients this implies with the asymptotic 



formula (2.2) that one gets in the first case just standard damped oscillations of har- 



monic type. In the latter case, several singularities in the complex plane approach the 
real axis which leads to more complicated oscillations (several phases) as can be seen 
in Fig. |2.17| Fitting the quantity in (p3|, we get a = 2.1692 and C = 241.7659. 



The additional structure in the oscillations of the Fourier coefficients implies that a 
linear fitting of the phase 4> is not efficient as can be seen in Fig. 2.17 



If we trace the quantities a and S during the computation as a function of time 
for k > 10 (again there are not enough non-trivial Fourier coefficients for k > 100 
at small times), we get the expected behavior as can be seen in Fig. 2.18 The used 
fitting just takes care of the singularity in the complex plane closest to the real axis. 
The imaginary part 5 of the singularity decreases strongly with t before the critical 
time of the Hopf solution and comes close to the axis there. After t c the decrease 
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Fig. 2.17. Difference between the mod ulus of the Fourier coefficients for the situation of 
Fig. 2.l6\ and the fitted asymptotic formula \2.2ty for k > 10 on the left, and for the quantity <f> 
12. bT) on the right. 



becomes very slow, and the singularity stays at a finite distance from the real axis. 
But again there is no clear indication of the precise values of the critical quantities. 
The real part of the singularity a can also be seen in Fig. |2.18| It gives in some sense 
the almost constant speed of the dispersive shock. 





Fig. 2.18. Fitting parameters for the solution to the KdV equation for the initial data uq(x) 
sech 2 x for e = 0.01 for S on the left, and a on the right. 



3. Kadomtsev-Petviashvili equations. In this section we will study shock 
formation in the dKP equation with the methods discussed for the example of the 
Hopf equation in the previous section. This will be done for both the dKP I and 
dKP II equations for two classes of initial data. The first class is motivated by the 
line solitons of KP, i.e., solutions exponentially localized in one spatial dimension 
and infinitely extended in the other. The second class of initial data are localized in 
both spatial directions. Then we will study for several values of the small dispersion 
parameter e how the difference between the dKP solution and the corresponding KP 
solutions scales at the critical point and for t -C t c in dependence of e. 

3.1. Theoretical preliminaries. We fist collect some known analytic facts on 



KP and dKP equations we will need in the following. The KP equation (1.5 1 is not 



in the standard form for a Cauchy problem since t is not a timelike coordinate if it 
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is discussed as a standard second order PDE, see e.g. [25 . To be able to treat initial 



value problems in t, equation (1.5) can be rewritten in evolutionary form, 



dtu + 6ud x u + e 2 d xxx u = —Xd x 1 d yy u, A = ±1. (3-1) 



J yy L 
5-1 



The antiderivative d x in (3.1 1 is to be understood as the Fourier multiplier with the 



singular symbol —i/k x . In the numerical computation this multiplier is regularized 



as described in [53]. Equations (3.1) and (1.5) are equivalent for certain classes of 



boundary conditions as periodic or rapidly decreasing at infinity, i.e., the ones studied 



in the present paper. The evolutionary form of dKP follows from (3.1 ) for e = 



The divergence structure of the KP equations (1.5) has the well known conse- 
quence that 



/ 

J 7 



d yy u(x, y, t)dx = 0, Vi > 0, (3.2) 



vv 

T 



even if this constraint is not verified at the initial time. As shown in [141 I33j the 
solution to a Cauchy problem not satisfying the constraint will not be smooth in time 
for t — 0. This leads to numerical problems which could have a negative impact on the 
determination of break-up singularities in dKP equations, see |25j . Thus we consider 
only initial data subject to the constraint (3.2). The first example with period 2irL y 



in y is 

u{x, y, 0) = exp (-(a: - cos(y/L y )) 2 ) , (3.3) 

whereas the second is just the derivative of a rapidly decreasing (in both spatial 
dimensions) function as studied in [25l [24] , 



u(x,y,0) = -d x scch 2 {R), R = y/x 2 +y 2 . (3.4) 

It is known (see e.g. [35] for references and examples) that solutions to KP equations 
for initial data in the Schwarzian space of rapidly decreasing functions generically do 
not stay in this space, but develop tails with an algebraic fall off to infinity. These 
tails will lead to a Gibbs phenomenon in a periodic setting. To reduce the latter, we 
will choose a larger computational domain than would be necessary if the solutions 
remained in the Schwarzian space as for instance in the case of KdV. We always give 
the Fourier coefficients to show that these algebraic tails do not affect the results in 
our examples. 

The dKP equation has only three conserved quantities and thus does not be- 
long to the family of completely integrable equations having an infinite number of 
conserved quantities. Therefore standard solution generating techniques as dressing 
transformations |34j and linear Riemann-Hilbert problems cannot be applied in this 
context. In [12] it was shown that solutions to the dKP equation can be found in 
terms of Einstein- Weyl geometries which can be constructed with Twistor methods, 
see for instance [35J H3]. However this approach is rather implicit if one wants to 
construct solutions for a given Cauchy problem and has not yet been used to this end. 
In [13] it was shown that the dKP equation allows for an infinite number of hydro- 
dynamic reductions and is thus integrable in this sense. Alinhac and coworkers [2] 
found that the dKP equation appears as a universal model in the geometric analysis 
of blow-up in nonlinear wave equations. They showed that finite time wave breaking 
is generic for solutions to this equation, proved existence of regular solutions up to the 
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break-up time t c and classified the singularity. Localized initial data become singular 
in a point (t c , x c , y c ) and the type of singularity is a one-dimensional cusp as for the 
Hopf equation, the second spatial dimension remains regular. However these methods 
do not provide formulae to determine (t c ,x c ,y c ) for given initial data. 

The most explicit results for break-up in dKP solutions so far have been given by 
Manakov and Santini [28] . They showed in [30] that the dKP equation can be solved 
in implicit form 

u = F(x-2ut,y,t), (3.5) 

similar to the solution of the Hopf equation via the method of characteristics. But F 
is here an integral over one component of the solution of a nonlinear vector Ricmann- 
Hilbcrt problem, see [28] for details. The latter is equivalent to a systems of two 
coupled nonlinear, singular integral equations. If the latter system of integral equa- 
tions is solved for a given normalization condition, the dKP solution follows from the 



implicit relation (3.5). To get a numerical solution to dKP in this way is not straight 
forward and has not been achieved for general initial data so far. But the form of the 
solution allows the study of singularity formation which has been done in [35] |2"§] . 
It is shown that there will be generically a gradient catastrophe at (t c ,x c ,y c ) in one 
spatial direction, whereas the solution remains smooth in the second spatial direction. 
Formulae are given for the solution for small \t — t c \, \x — x c \ and \y — y c \, and the type 
of singularity is identified again as cubic in the singular direction. In particular, they 
show that at the breaking time t c , if the initial condition uq(x, y) for the dispersionless 
KP equation is even in y, then the first breaking point is on the x-axis, i.e., y c = 0, 
and that the y-derivative of u is zero at (x c , 0). Our numerical results for the case of 
localized initial data in Sect. 3.3. are in accordance with this prediction. But it is not 
clear how to obtain in the formalism [25] the values for the critical point (t c ,x c ,y c ) 
for given initial data. Therefore we solve the dKP equation in this paper directly. 

The above mentioned theoretical results in [2] [28] and references therein indicate 
that break-up in dKP solutions is one-dimensional. We choose initial data symmetric 
with respect to the y-coordinate to ensure a blow-up of the x-gradient only. However 
this is without loss of generality since the same approach as applied below can be 
used for general initial data. In this case one first has to determine numerically the 
direction of the strongest gradient which will indicate where break-up occurs. Since 
the rest is as for the symmetric initial data, we concentrate on this conceptionally 
simpler case. 

3.2. Numerical approach. We use here again a Fourier spectral method for the 
spatial coordinates and as for the Burgers' and KdV equations in the previous section 
a fourth order exponential time differencing scheme derived by Cox and Matthews in 
[7] for the time integration. We consider periodic (up to numerical precision) solutions 
in x and y, i.e. solutions on T 2 x R. The computations are carried out with N x x N y 
points for (x,y) e [—L x w, L x ir] x [—L v tt, L v tt]. 

As discussed in the previous section, for a reliable identification of the shock 
formation via the asymptotic behavior of the Fourier coefficients, sufficient spatial 
resolution is crucial. In practice we needed 2 14 to 2 16 modes for the Hopf equation in 
1 + 1 dimensions. To allow such high resolution simulations, we had to parellelize the 
codes also for the asymptotic analysis of the Fourier coefficients. A prerequisite for 
parallel numerical algorithms is that sufficient independent computations can be iden- 
tified for each processor, that require only small amounts of data to be communicated 
between independent computations. To this end, we perform a data decomposition, 
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which makes it possible to do basic operations on each object in the data domain 
(vector, matrix...) to be executed safely in parallel by the available processors. 

Our domain decomposition is implemented by developing a code describing the 
local computations and local data structures for a single process. Global arrays are 
divided in the following way: denoting by x n = 2imL x /N x , y m = 2TrmLy/N y , n = 
—N x /2,...,N x /2, m — —Ny/2,...,N y /2, the respective discretizations of x and y in 
the corresponding computational domain, u is then represented by an N x x N y matrix. 
The latter is distributed among processors such that each processor Pi,i — l...n p , (n p 
denoting the number of processors we can access) will receive N x x — - elements of u 
corresponding to the elements 

/v. /v. \ 
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in the global array, and then each parallel task works on a portion of the data. 

While processors execute an operation, they may need values from other proces- 
sors. The above domain decomposition has been chosen such that the distribution of 
operations is balanced and that the communication is minimized. The access to re- 
mote elements has been implemented via explicit communications, using sub-routines 
of the MPI library. The computation of the 2d-FFT is performed by using a transpo- 
sition approach which takes advantage of the existing sequential FFT routines (more 
precisely, we use serial FFTW routines 16J) . The asymptotic fitting of the Fourier 
coefficients, here in one spatial direction, requires in addition 2 local communications, 
a core being reserved to do the required computations. More precisely, the proces- 
sor which 'takes care' of v(k x ,0,t) sends these data to the reserved processor, which 
receives the data (thus 2 local communications per time iteration, a SEND, and a 
RECEIVE), and then it alone performs the least square fit required, which takes less 
time than the computation of the solution for one iteration, performed by the other 
processors. The main technical problem in the use of ETD schemes is the efficient 
and accurate numerical evaluation of the functions 

<j>i(z) = -^— / e^V^dr, i = 1, 2, 3, 4, 

i.e., functions of the form (e z — l)/z and higher order generalizations thereof, as 
explained in |22j and in |39j . We use the approach given in [39) in our implementation. 
The computation of these functions takes only negligible time for the 2+1-dimensional 
equations studied here, especially since it has to be done only once during the time 
evolution. In addition each processor only computes the portion of these functions 
that the further local computations require, according to the domain decomposition 
specified above. 

To control the numerical accuracy we trace as for the Hopf equation in the pre- 
vious section and as in |24j the conservation of the L^ norm of u, the mass of the 
solution. It is exactly conserved for the KP equations, but due to unavoidable nu- 
merical errors the computed mass m(t) will numerically depend on time. Since mass 
conservation is not implemented in the code, it provides a valid check of the numer- 
ical accuracy for sufficient resolution in Fourier space. It was shown in [24| that the 
quantity 

m(t) - m(0) 

171(0) 
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typically overestimates the L^ norm of the difference between numerical and exact 
solution by two orders of magnitude. 

3.3. Shock formation in dKP solutions for initial data localized in one 
spatial dimension. We will first study numerically the appearence of break-up sin- 



gularities in solutions to the dKP equations for the initial data (3.3) with the methods 
explained above. Whereas initial data localized in both spatial dimensions will have 
a gradient catastrophe in just one point, initial data localized in one dimension and 
infinitely extended in the other appear to blow up on a curve. 

use N x = N y = 2 14 , L x = 

u 



To determine a solution for the initial data (3.3|, we use N x — A 
j y = 5. We observe that the solution develops a shock in the x-direction at time 
t ~ 0.19, as can be seen in Fig. |3.1[ 
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Fig. 3.1. Solution to the dKP I equation for the initial data 13. 31), for several values of t 



To determine the time of the appearance of this shock, we apply the fitting 
procedure described in the previous section for the Fourier coefficients (again de- 
noted by v{k x ,ky,t)) in ^-direction, i.e., v(k x ,0,t), to formula (2.3). As before (see 



Section 2.4) we use at least half of the Fourier coefficients with values above the 
rounding error, and an appropriately chosen interval [k m i n ,k max \. In this way, for 
N x = 2 14 we find the minimal possible precision to be p = 0.01, for k m i n — 5 and 
k-max — max(fc :E )/2. For this prescribed error, we find that t c ~ 0.1934, where 5 
vanishes, and B{t c ) = 1.337 ~ |. We observe at this time the formation of a shock 
in the ^-direction, which can be seen in Fig. 3.2 where we show the behavior of the 
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solution at different times, plotted on the x-axis, and the corresponding Fourier co- 
efficients, plotted on the fc^-axis. The time dependence of the fitting parameters is 





* t=0.05 

* t=0.10 

* t=0.15 

. t =0.1934 



1000 1500 



FlG. 3.2. Solution to the dKP I equation for the initial data \3.3\j , plotted on the x-axis for 
several values oft, and the corresponding Fourier coefficients, plotted on the k x -axis. 



shown in Fig. 3.3 As expected, we observed a rapid decrease of S(t), and the change 











2 




• — m 
6=0 


- 


1-5 


I 


0.5 


^^^ 










0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 



FlG. 3.3. Time dependence of the fittin g parameters 8 (left), B (right) for the numerical solution 
to the dKP I equation for initial data \3. 3U . The fitting is done for 5 < k x < max(fc I )/2. 



of B(t) from a value ~ 1.5 to the value ~ 1.33 when approaching the critical time, 



determined by the vanishing of 8(t). In Fig. 3.4 we show the time dependence of the 



fitting parameters close to t c , for different values of p. By doing a fitting of Slog(u), 
we find that at t c , x c = a(t c ) = 2.4134. 

In this case, we reach an error in the fitting of p ~ 0.01, which is the same as 
obtained for the Hopf equation. Note also that the approximation of B is of the same 
precision as in the latter. However, such a precision cannot always be achieved as will 
be seen below. 



The time evolution of the quantity A^ in (3.7), used as an indicator of the 
numerical accuracy, and of the Loo-norm of d x u until the critical time are shown 
in Fig. 3.5 The quantity A^ indicating the numerical mass conservation increases 
close to t c , but stays below 10~ 7 , which implies that the system is still well resolved 
up to the shock formation. Whereas the L^ norm of u stays finite, there are strong 
indications of a blow-up of u x close to the critical time, which indicates that the fitting 
in fact identified correctly the gradient catastrophe. 
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Fig. 3.4. Time dependence of the fitting parameters 8 (left), and B (right) close to t c ~ 0.1934 
for the numerical solution to the dKP I equation for the initial dat a \ 3. 3\ for different values of p. 
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Fig. 3.5. Time evolution of the quantity Ag as an indicator of the numerical accuracy (left), 
and of Hc^^Hoo (right) for the solution to the dKP I equation for initial data 



As one can see in Fig. 3.6 where \d x u\ and \d v u\ are plotted at t = 0.1934, 



the gradient catastrophe does not appear in only one spatial point for the infinitely 
extended initial data. 





Fig. 3.6. Derivatives of the dKP I solution for the initial data \3.3\ at the critical time t c 
0.1934, \d x u\ (left) and of\d y u\. 



The situation is rather similar for the dKP II case, for the same initial data (3.3) 



By using the same fitting bounds as before, (k m i n = 5 and k max = max(fc)/2), we 
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find that the solution of the dKP II equation develops a shock at t c ~ 0.1934, which 
is exactly the previously identified value for the critical time for dKP I. The precision 
we could achieve here is also similar, A = 0.01, B(t c ) = 1.337 and a(t c ) = 2.41. The 
solution at t c — 0.1934 is shown in Fig. 3.7 







Fig. 3.7. Solution to the dKP II equation for the initial data (S.Slf, at t c = 0.1934 



A noticeable difference appears however in the profile of u x . One can observe in 
Fig. |3.8[ that in contrast to the dKP I case, the strongest value of the cc-gradient here 
occur at the trailing part of the solution, whereas in the case of dKP I, it was occurring 
at roughly y = 0. This is related to the focusing and defocusing character of dKP I 
and dKP II respectively. In the first case, the most advanced point in propagation 
direction gets focused, whereas in the latter case, this happens for the trailing points. 





Fig. 3.8. Derivatives of the dKP II solution for the initial data 1373V at the critical time 
t c = 0.1934, \8 x u\ (left) and of \8 y u\. 



3.4. Shock formation in dKP solutions for initial data localized in both 
spatial dimensions. In this subsection we will study shock formation in dKP solu- 



tions for initial data of the form (3.4) which are localized in both spatial directions. 
It is known, see [5J [2S], that localized data with a single maximum will lead to a 
break-up in a single point as expected. 

Remark 3.1. In order to satisfy the constraint (3.2), the data studied here do 
not have a single maximum. Therefore there will be actually two points of gradient 
catastrophe. The second gradient catastrophe appears for dKP I for negative x for a 
slightly larger t (for dKP II the role of the two critical points is interchanged) . Since 
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the code can only be reliably run until the first shock formation, we cannot obtain the 
second break-up with the method, and we will only discuss this time in the following. 

= 5, we observe that the solution develops 
0.22, as we can see in Fig. 



For N x = N y = 2 i4 and L x = L, 
a shock in the x-direction at a time t 
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Fig. 3.9. Solution to the dKP I equation for initial data (3. 4\), for several values of t 



accordance with the experiments in |25| for the same initial data, where the break-up 
was studied qualitatively for a dissipative regularization of the dKP equation. There 
it was found that the gradient catastrophe for x > is reached roughly at t ~ 0.23. 
In Fig. |3.9| it can also be seen that the initially rapidly decreasing function u develops 
tails with an algebraic fall-off to infinity. Due to the imposed periodicity condition, 
this leads to a Gibbs phenomenon at the boundary of the computational domain. 
To reduce the effect of this phenomenon on the Fourier coefficients, we have used 
a considerably larger computational domain than shown in Fig. |3.9| We present in 
Fig. 3.10 the solution to the dKP I equation with initial data (3.4), plotted on the 



x-axis for several values of t and the corresponding Fourier coefficients, plotted on 
the fc^-axis. From the latter it can be inferred that the asymptotic behavior of the 
Fourier coefficients is dominated by the break-up singularity and not by the algebraic 
fall off of the solution. 

In this case, we find that the minimal attainable fitting error (with the above 
described requirement to use at least half of the Fourier coefficients) is p ~ 0.5 for 
k m in = 10 and k max = max(k x )/2. For such bounds, 6 vanishes at t c = 0.2216, 
where B(t c ) = 1.34 and A = || hx\v{k x ,Q)\ - (A - B\nk x - kj)^ = 0.48. The 
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Fig. 3.10. Solution to the dKP I equation for initial data \3.4\ , plotted on the x-axis for several 
values oft, and the corresponding Fourier coefficients, plotted on the k x -axis. 



considerably larger error in the fitting than for the initial data (3.3) is related to 



the appearence of a second break-up singularity for the data (3.4 1 as explained in 
Remark |3.1| Nonetheless the results of the fitting are reliable as can be seen in 
the following. A fitting of Qlog(i>) yields x c — a(t c ) — 1.7961. The value of B is 
compatible with the expected value | . The time dependence of the fitting parameters 



is shown in Fig. 3.11 for t ~ t c . 




Fig. 3.11. Time dependence of the fitting parameters 6 (left), B (right) for the numerical 
solution to the dKP I equation for initial data\3.4\ The fitting is done for 10 < k x < ma,x(k x )/2. 



The time dependence of the quantity A^ (3.7) controlling the numerical error, 
and of ||w X ||c 



until t c = 0.2216 are shown in Fig. 3.12 



In this case, the gradient catastrophe appears in only one point (x c , y c ), as we can 
see from Fig. 



3.13 



where we show the behavior of \d x u\ (left) and of \d y u\ (right) at 
t = 0.2216. Note that the gradient in y-direction is much smaller than in ^-directions 
which confirms the theoretical expectation that break-up only occurs in one direction 
whereas the solution remains regular in the other direction. This is even more obvious 



in a close up in Fig. 3.14 (top), where we show \d x u\ in the region where the gradient 
catastrophe occurs. We recover the previously computed value a(t c ) ~ 1.79 in Fig. 

1.79,0), and 
4, 



3.14 (bottom), where \d x u(x,0)\ is shown. It can be seen that (x c ,y c ) 

13 



we find that u y (t c ,x c , 0) = 4.8* 10 
in [29 for symmetry reasons. 



0, as also predicted by the theoretical results 



The situation is similar for the dKP II case for the initial data (3.4), which we 
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Fig. 3.12. Time evolution of the quantity Ag (3. 7| ) (left) and of \\u x \\oo (right) for the solution 
to the dKPI equation for the initial data 13. J^. 





Fig. 3.13. Der ivati ves \d x u\ and of \d y u\ at t = 0.2216 for the solution to the dKP I equation 
for the initial data \3. 4\ . 
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Fig. 3.14. Close up of Fig. \3.13\ of \d x u\ (left) and the profile of \d x u(x, 0)| (right) at t = 0.2216. 



treat analogously. We find that the fitting used for dKP I is optimal also for dKP 
II (k m i n = 10,k max = max(fc)/2), which yields the vanishing of 6 at t c ~ 0.2216, 
B(t c ) = 1.34, and A = || In \v(k x , 0)| - (A- B\nk x ~ k x 8)\\oo = 0.48. We observ e that 
the solution develops a shock in the x-direction (x < 0) at time t c , see Fig. 3.15 where 
we show the solution to the dKP II equation with initial data (3.4 1, at t = t c . Note 



that the tails with the algebraic fall off are now directed towards — oo in contrast to 
the dKP I case. Thus the gradient catastrophe appears in this case first for negative 
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x (again we only treat the point of gradient catastrophe appearing first). By doing a 
fitting of 51og(v), we find a(t c ) ~ —1.79. The good agreement of this value with the 
maximum of the gradient of u shows the self consistency of the approach. 





Fig. 3.15. Solution to the dKP II equation with initial data \3.4\ at t=0.2216 



In this case, the gradient catastrophe appears again only in one point (x c ,y c ) 



( — 1.79,0), as we can see from Fig. 3.16 where we show \d x u\ at t = 0.2216, and 
u y (t c ,x c ,0) = 3.3 * 10 -12 . The solutions to both dKP equations show thus a similar 



100 




Fig. 3.16. Derivative \d x u\ at t = 0.2216 of the solution in Fig. \3.15\ 



behavior, independent of the sign of A in (3.1 1, except that of course, in the case 
of dKP I, the gradient catastrophe appears at (x c ,y c ) = (1.79,0), whereas in the 
case of dKP II, it appears at (x c ,y c ) = (—1.79,0) (as actually expected). As already 
mentioned, this is due to the tails with an algebraic fall off towards infinity which are 
for dKP I directed towards +oo, whereas they are directed towards — oo for dKP II. 

3.5. KP solution in the small dispersion limit for initial data localized 
in one spatial dimension. Solutions to the KP equation will develop dispersive 
shocks, zones of rapid modulated oscillations, in the vicinity of a shock of solutions to 
the dKP equation for the same initial data. These were studied numerically for the 
first time in [35] , see also [21] • For KdV it was found that the difference between Hopf 

30 



and KdV solution scales as e 2 for t -C t c and as e 2 ' 7 for t ~ t c . In this subsection, 
we establish numerically such scaling laws for the initial data localized in only one 
spatial direction. 



For initial data of the form ( 3.3 ) , we show the solution to the KP I equation in the 



small dispersion limit with e = 0.1 in Fig. |3.17| for several times. The computations 
are carried out with 2 14 x 2 14 points for x x y E [— 5tt, 5ir] and A t — 4 * 10 -5 . 
As expected, the dispersive regularization of the break-up singularity leads to rapid 



t=0.1 
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Fig. 3.17. Solution of the KP I 
with e = 0.1 for several values oft 



equation in small dispersion limit for the initial data 1373\ 



modulated oscillations in the region where the corresponding dKP I solution develops 
a shock. 

The number of oscillations increases as e tends to 0, as already observed in |25j . 
This can be seen in Fig. |3.18| where the solution to the KP I equation in the small 
dispersion limit is shown for e = 0.02 at t — 0.4. In Fig. |3. 19] we present the solutions 
to the KP I equation in the small dispersion limit for several values of e on the x-axis 
for t = 0.4. 



The contourplots of the solutions at t — 0.4 are shown in Fig. 3.20 for several values 
of e. 

As usual, we ensure that the system is numerically well resolved by checking the 
decay of the Fourier coefficients, and the conservation of the numerically computed 
mass. The Fourier coefficients of the solution of the KP I equation in the small 
dispersion limit are plotted on the fc^-axis for several values of e in Fig. |3.21| at 
t = 0.4. They decrease to machine precision in all cases. The quantity A# is always 
smaller than 10~ 8 , which indicates a numerical error well below plotting accuracy. 

An important question is the scaling with e of the L^, norm of the difference 
between dKP and KP solutions for the same initial data. The L^ norm of this 

0.1 - t c /2 (left) and at t c = 0.1934 (right) 
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difference is shown in Fig. 3.22 



at t 




Fig. 3.18. Solution of the KP I equation in the small dispersion limit for initial data of the 
form JO) with e = 0.02 at t=04- 
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Fig. 3.19. Solutions to the KP I equation in the small dispersion limit on the x-axis att = 0.4 
for several values of e 



in dependence of e for 0.01 < e < 0.1. A linear regression analysis (log 10 Aoo = 
alog 10 e + b ) shows that A^ decreases as 

O (e 196 ) - O (e 2 ) at t = 0.1 ~ i c /2, with a = 1.961 and b = 0.6840 (3.8) 
(e 028 ) ~ O (e 2/7 ) at t = t c = 0.1934, with a = 0.278 and b = -0.4214. (3.9) 

In both cases, the correlation coefficient is r = 0.999. 

The behavior of solutions of KP II in the small dispersion limit, for initial data of 
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Fig. 3.20. Contour plots of the solutions to the KP I equation in the small dispersion limit at 
t = 0.4 for several values of e. 





Fig. 3.21. Fourier coefficients to the solution of the KP I equation in the small dispersion limit 
on the k x -axis for several values of e, at t = 0.4 (left), and the time evolution of the quantity Ag 
fO) (right) 



the form (3.3 1 is completely similar to the KP I case, see for example Fig. 3.23 where 



we show them plotted on the s-axis, for different values of e. Note that the oscillations 
are somewhat smaller than in the KP I case because of the discussed defocusing effect 
of KP II. We repeat the same study as before, and find that also in this case, the 
system is well resolved until t = 0.4, and that the L^ norm of the difference between 



dKP II and KP II solutions for initial data ( 3.3 ) decreases as in the case of dKP I/KP 
I: We find that Aqq decreases as 

O (e 1 ' 96 ) - O (e 2 ) at t = 0.1 ~ t c /2, with a = 1.9608 and b = 0.6840 (3.10) 

O (e 028 ) - O (e 2,7 \ at t = t c = 0.1934, with a = 0.2751 and b = -0.4122. (3.11) 

The correlation coefficient is r = 0.999 for both times. 

33 




Fig. 3.22. Loo norm Aoo of the difference between KP and dKP solution for the initial data 
3.3ty in dependence of e at t = 0.1 ~ t c /2 (left) and at t c = 0.1934 (right) for several values of e. 
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Fig. 3.23. Solutions to the KP II equation in the small dispersion limit on the x-axis att = 0.4 
for several values of e 



3.6. KP solution in the small dispersion limit for initial data localized 
in both spatial dimension. The same study as in the previous subsection will be 



carried our for the initial data (3.4). Interestingly the results are very similar though 



the (first) dKP shock appears here in just one point. 



For the initial data of the form (3.4) one gets the dispersive shock already studied 
in [55] ■ The solution at time t = 0.4 > t c can be seen in Fig. 3.24 for several values 



of e. It can be recognized that the number of oscillations increases with decreasing e, 
and that the oscillations are more and more confined to a well defined zone. This is 
even more visible on the x-axis as can be inferred from Fig. |3.25| The corresponding 
contour plots are shown in Fig. |3.26| 



The Fourier coefficients of the solution on the fc^-axis can be seen in Fig. 3.27 for 
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Fig. 3.24. Solutions to the KP I equation in the small dispersion limit for t = 0.4 for several 
values of e 



several values of e at t = 0.4. They decrease to machine precision in all cases. The 
quantity A# used as an indicator for the numerical accuracy is always smaller than 
10 -9 , see Fig. 3.27 which is more than satisfactory for our purposes. 

For the initial data (3.4), we find a similar scaling with e as for the data (3.3): We 
obtain via a linear regression analysis (log 10 Aoo = a log 10 e + b ) that the L M norm 
Aqo of the difference between dKP and KP solution decreases as 

O (e 194 ) - O (e 2 ) at t = 0.125 < t c , with a = 1.943 and b = 1.2239 (3.12) 
O (e 030 ) ~ O (e 2/7 \ at t = t c = 0.2216, with a = 0.30 and b = -0.5122. (3.13) 

In both cases, the correlation coefficient is r = 0.999. 

The results for the KP II case are very similar to the KP I. We obtain again 
a dispersive shock in the small dispersion limit, see Fig. |3.30| and Fig. |3.29[ already 
observed in [2S] for the same initial data. Since the tails with the algebraic fall off are 
now directed towards — oo, the first shock for dKP as well as the stronger oscillation 
appear for negative x. The mass transfer towards — oo simply implies that there 
is more mass in this region which triggers a dispersive shock. The corresponding 
solutions on the x-axis can be seen in Fig. |3.30| The related contour plots are shown 
in Fig. |3.31| The Fourier coefficients of the solution to the KP II equation in the 



small dispersion limit are plotted on the k x -axis for several values of e in Fig. 3.32 
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Fig. 3.25. Solutions to the KP I equation in the small dispersion limit on the x-axis for t = 0.4 
for several values of e 
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Fig. 3.26. Contour plots of the solutions to the KP I equation in the small dispersion limit for 
t = 0.4 for several values of e 



at t — 0.4. They decrease to machine precision in all cases, which ensures that the 
system is well resolved until t = 0.4. The conservation of the numerically computed 
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Fig. 3.27. Fourier coefficients to the solution of the KP I equation in the small dispersion limit 
on the k x -axis for several values of e, at t = 0.4 (left) and the time evolution of the quantity Ag 
fO) (right) 




Fig. 3.28. Error A c 
for several values of e. 



in dependence of e at t = 0.125 ~ t c /2 (left) and at t c = 0.2216 (right) 



mass always reaches a precision of better than 10 -9 , see Fig. 3.32 

The scaling laws are also almost identical to the previous case: We find, with a 
linear regression analysis (log 10 A^ = a log 10 e + b) that A^ decreases as 

O (e 198 ) - O (e 2 ) at t = O.f 1 - i c /2, with a = 1.9815 and b = 0.8998 (3.14) 
O (e 029 ) - O (e 2/J ) at t = t c = 0.2216, with a = 0.295 and b = -0.5288. (3.15) 

In both cases, the correlation coefficient is r — 0.999. 

Thus in all cases, we find approximately the same scaling laws as in the KdV 
case. 

4. Outlook. In this paper we have shown that it is possible to use asymptotic 
Fourier analysis to determine the time of the singularity formation in nonlinear evo- 
lution equations. A precondition for the applicability of the method is obviously that 
the PDE is solved with sufficient accuracy and resolution. The former can be con- 
trolled via a conserved quantity, typically the L2 norm or the energy, the latter via the 
decrease of the Fourier coefficients for large wave numbers. The order of magnitude 
of the Fourier coefficients for the highest wave numbers is in some sense a measure for 
the numerical resolution. The conserved quantities, which will numerically depend 

37 



£=0.08 



£=0.06 





1 , 




Fig. 3.29. Solutions to the KP II equation in the small dispersion limit at the maximal time 
of computation, t = 0.4 for several values of e 
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Fig. 3.30. Solutions to the KP II equation in the small dispersion limit plotted on the x-axis 
at the maximal time of computation, t = 0.4 for several values of e 



on time because of unavoidable errors, cannot indicate a higher precision than this 
resolution, and typically overestimate the accuracy by 2-3 orders of magnitude. If 
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Fig. 3.31. Contour of the solutions to the KP II equation in the small dispersion limit at the 
maximal time of computation, t = 0.4 for several values of e 
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Fig. 3.32. Fourier coefficients to the solution of the KP II equation in the small 
limit, plotted on the k x -axis for several values of e, at t = 0.4 (left) and time evolution of the mass 
conservation Ag (right) 



sufficient numerical resolution is provided, normally more than 2 14 Fourier modes, a 
minimal wave number k m i n for the least square fitting to the asymptotic formula (2.2 1 



has to be chosen. This has to be clearly larger than 1, but should include sufficient 
Fourier coefficients with larger modulus than the rounding error for the considered 
times. Values between 10 and 100 appear to be convenient here. Then one concludes 
an aimed at fitting error in dependence of the resolution (several runs might be nec- 
essary to make this self consistent, but generally it can be chosen of the same order 
as the estimated accuracy of the solution) . The upper limit for the wave numbers for 
which the fitting is performed is then chosen as the maximal |fe| for which the dif- 
ference between fitted curve and Fourier coefficients on the whole considered interval 



:-!<) 



is smaller than the numerical accuracy of the solution. If this cannot be achieved 
with at least half of the Fourier coefficients with a modulus larger than the rounding 
error, the number of Fourier modes has to be increased. In this way it is possible to 
trace the singularity with acceptable precision up to the point where it hits the axis, 
which defines the critical time. This approach gives the critical quantities with the 
necessary precision to study the scaling of solutions to a dispersive regularization as 
KdV for Hopf. 

In 2 + 1 dimensions and higher, the needed resolution can only be achieved on 
parallel computers. As an example we studied shock formation in the dKP equation 
which could be treated with a parallel version of the code with the same precision 
as for Hopf for initial data which are infinitely extended in one direction. For data 
being localized in all spatial directions, the achievable accuracy was slightly lower. 
The reason for this is the appearance of a second singularity shortly after the first 
we studied. This is due to the fact that the initial data for KP have to satisfy the 



constraint (3.2) in order to allow a smooth solution in time for t = 0. The presence of 
a second singularity in the complex plane limits somewhat the accuracy of the fitting 
procedure for the Fourier coefficients. Nonetheless we could identify the critical point 
and the critical singularity with the needed precision. This allowed to study the 
scaling of the difference between KP solutions in the small dispersion limit and the 
dKP solution at break-up. It was shown that the same scaling (within numerical 
tolerance) is observed at the critical time as in the KdV case, e 2 / 7 . The task is now to 
give an asymptotic description of the KP solutions in the small dispersion limit close 
to the break-up of dKP as in [5] . 

The numerical analysis presented in this paper for dKP will also be applied to 
other 2+1- and higherdimensional PDEs, for instance the Davey-Stewartson equations 
in the semiclassical limit. One problem in this context is that the singularity might 
in these cases be genuinely multi-dimensional, not as for dKP where the gradient 
blows up only in one spatial direction. The essentially one-dimensional approach can 
of course still be applied in this case. Alternatively instead of the Fourier transform 
u(k,t), one can consider as in |42j the angle averaged energy spectrum defined by 

E{K,t)= ]T \u(k',t)\ 2 , 

K<\k'\<K+l 

where \k'\ — wfc^ + fc 2 Slightly weaker estimates hold for E(K,t) for an analytic 
function u in S a , and thus, to apply the asymptotic fitting to the Fourier coefficients, 
one assumes that E{K, t) — C (t)K~ a ^ e~ s ^ K . Similary to the one-dimensional case, 
the appearance of a real singularity implies that 5(t) vanishes at a finite time t c . This 
implies that sufficient numerical resolution has to be provided which can as in the 
present paper in general only be achieved on parallel computers. With this approach 
it should be possible to study singularity formation in general nonlinear evolution 
equations without dissipation and dispersion. 
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